!#############################################################
 module markeraux
! Auxliliary routines used in combination with the markers 
! or for interpolation
!#############################################################

!*******************
  contains
!*******************
 

!=====================================================================================
 subroutine grid2mk(xm,zm,grdval,mkval)
! Given the grid, find the four nodes that surround (xm,zm) 
! and compute the value 'mkval' via linear interpolation from grid values 'grdval' 
!
! Input: 
! xm => real. x-coord of the interpolation point
! zm => real. z-coord of the interpolation point
! grdval => real, array(0:nx+1,0:nz+1). Field to be interpolated
! Output:
! mkval => real. Interpolated value at (xm,zm)
!=====================================================================================
   use share; use domain
   implicit none
   integer :: i1,i2,j1,j2
   real (kind=8), intent(in) :: xm,zm
   real (kind=8), intent(out) :: mkval
   real (kind=8) :: dxm,dzm
   real (kind=8), dimension(0:,0:) :: grdval
!===========================================================
 
!-----------------------------------------
! i1,j1 => integer. lower left grid point
! i2,j1 => integer. lower right grid point
! i1,j2 => integer. upper left grid point
! i2,j2 => integer. upper right grid point
!-----------------------------------------
   i1=bisec(xgrid,xm)
   i2=i1+1
   dxm=xm-xgrid(i1)   
   if(dxm.eq.0d0) i2=1
   
   j1=bisec(zgrid,zm)
   j2=j1+1
   dzm=zm-zgrid(j1)
   if(dzm.eq.0d0) j2=1
   
! Interpolate the field 'grdval' to the marker
   mkval=(1-dxm/delx(i2))*(1-dzm/delz(j2))*grdval(i1,j1) &
        +(dxm/delx(i2))*(1-dzm/delz(j2))*grdval(i2,j1)   &
        +(1-dxm/delx(i2))*(dzm/delz(j2))*grdval(i1,j2)   &
        +(dxm/delx(i2))*(dzm/delz(j2))*grdval(i2,j2)

 end subroutine grid2mk
 

!===================================================================
 subroutine Wmarker(xm,zm,i,j,w)
! Find the four nodes that surround a marker (xm,zm) 
! and compute marker's statistical weight w
!===================================================================
   use share; use domain
   implicit none
   integer :: i1,i2,j1,j2
   integer, intent(in) :: i,j
   real (kind=8), intent(in) :: xm,zm
   real (kind=8), intent(out) :: w
   real (kind=8) :: dxm,dzm
!===================================================================
      
   i1=bisec(xgrid,xm)
   i2=i1+1
   if(xm.le.xgrid(i)) then
     dxm=xgrid(i2)-xm   
   else
     dxm=xm-xgrid(i1)
   endif
   
   j1=bisec(zgrid,zm)
   j2=j1+1
   if(zm.le.zgrid(j)) then
     dzm=zgrid(j2)-zm
   else
     dzm=zm-zgrid(j1)
   endif
      
! Compute the statistical weight for the marker
   w=(1d0-dxm/delx(i2))*(1d0-dzm/delz(j2))/delx(i2)/delz(j2)

 end subroutine Wmarker      
 

!==========================================================================
 subroutine Wmarker4(xm,zm,ai,aj,w)
! Compute marker's statistical weighs
!
! Given the position (xm,zm) of a marker, the routine finds the 
! four nodes that surround the marker and computes the marker's
! statistical weight with respect to each node
!
! Input:
! xm => real. Marker x-coordinate
! zm => real. Marker z-coordinate
! Output:
! ai => integer, array(4). Array of x indices
! aj => integer, array(4). Array of z indices
! w => real, array(4). Array of marker's weights
!==========================================================================
   use share; use domain
   implicit none
   integer :: i1,i2,j1,j2
   integer, dimension(4), intent(out) :: ai,aj
   real (kind=8), intent(in) :: xm,zm
   real (kind=8), dimension(4), intent(out) :: w
   real (kind=8) :: dxm1,dxm2,dzm1,dzm2,d
   real (kind=8), parameter :: dmin=1d-12
!===================================================================
      
!---------------------------------------------------------------
! Coordinate indices of the four nodes that surround the marker
! i1,j1: lower left grid point
! i2,j1: lower right grid point
! i1,j2: upper left grid point
! i2,j2: upper right grid point
!---------------------------------------------------------------
   i1=bisec(xgrid,xm)
   i2=i1+1
   dxm1=xm-xgrid(i1)
   dxm2=xgrid(i2)-xm   
   ai(1)=i1
   ai(2)=i2
   ai(3)=i2
   ai(4)=i1
   
   j1=bisec(zgrid,zm)
   j2=j1+1
   dzm1=zm-zgrid(j1)
   dzm2=zgrid(j2)-zm
   aj(1)=j1
   aj(2)=j1
   aj(3)=j2
   aj(4)=j2

!------------------------------------
! Compute markers wheights
!------------------------------------
   if(mkinterp.eq.1) then

!   Inverse distance
     d = dsqrt(dxm1**2d0 + dzm1**2d0)
     if(d.lt.dmin) d = dmin
     w(1) = 1d0/d

     d = dsqrt(dxm2**2d0 + dzm1**2d0)
     if(d.lt.dmin) d = dmin
     w(2) = 1d0/d

     d = dsqrt(dxm2**2d0 + dzm2**2d0)
     if(d.lt.dmin) d = dmin
     w(3) = 1d0/d

     d = dsqrt(dxm1**2d0 + dzm2**2d0)
     if(d.lt.dmin) d = dmin
     w(4) = 1d0/d

   elseif(mkinterp.eq.2) then

!   Statistical weight
     w(1) = (1d0-dxm1/delx(i2))*(1d0-dzm1/delz(j2))
     w(2) = (1d0-dxm2/delx(i2))*(1d0-dzm1/delz(j2))
     w(3) = (1d0-dxm2/delx(i2))*(1d0-dzm2/delz(j2))
     w(4) = (1d0-dxm1/delx(i2))*(1d0-dzm2/delz(j2))

   endif

 end subroutine Wmarker4     
 

!===================================================================
 subroutine rk2(xold,zold,xnew,znew)
! 2nd-order in space/1st order in time Runge-Kutta to
! find the new position (xnew,znew) of a marker initially 
! at (xold,zold) passively advected by the flow (u,v)
!===================================================================
   use share; use domain
   implicit none
   real (kind=8), intent(in) :: xold,zold
   real (kind=8), intent(out) :: xnew,znew
   real (kind=8) :: x1,x2,z1,z2,k1,k2
!===================================================================
   
   x1=xold; z1=zold
   call grid2mk(x1,z1,u,k1)
   x2=xold+dt/2*k1
   call grid2mk(x2,z1,u,k2)
   xnew=xold+dt*k2

   x1=xold; z1=zold
   call grid2mk(x1,z1,v,k1)
   z2=zold+dt/2*k1
   call grid2mk(x1,z2,v,k2)
   znew=zold+dt*k2

 end subroutine rk2   
 
!===================================================================
 subroutine rk4(xold,zold,xnew,znew)
! 4th-order in space/1st order in time Runge-Kutta to
! find the new position (xnew,znew) of a marker initially 
! at (xold,zold) passively advected by the flow (u,v)
! If one of the four computation points needed lies outside 
! the domain, the new position is found by a 2nd order Runge-Kutta
!===================================================================
   use share; use domain
   implicit none
   real (kind=8), intent(in) :: xold,zold
   real (kind=8), intent(out) :: xnew,znew
   real (kind=8) :: x1,x2,x3,x4,z1,z2,z3,z4,k1,k2,k3,k4
!===================================================================
   
   x1=xold; z1=zold
   call grid2mk(x1,z1,u,k1)
   x2=xold+dt/2*k1
   call grid2mk(x2,z1,u,k2)
   x3=xold+dt/2*k2
   if(x3.gt.dimx.or.x3.lt.0d0) then
     xnew=xold+dt*k2
   else
     call grid2mk(x3,z1,u,k3)
     x4=xold+dt*k3
     if(x4.gt.dimx.or.x4.lt.0d0) then
       xnew=xold+dt*k2
     else
       call grid2mk(x4,z1,u,k4)
       xnew=xold+(k1+2*k2+2*k3+k4)*dt/6
     endif  
   endif

   x1=xold; z1=zold
   call grid2mk(x1,z1,v,k1)
   z2=zold+dt/2*k1
   call grid2mk(x1,z2,v,k2)
   z3=zold+dt/2*k2
   if(z3.gt.dimz.or.z3.lt.0d0) then
     znew=zold+dt*k2
   else
     call grid2mk(x1,z3,v,k3)
     z4=zold+dt*k3
     if(z4.gt.dimz.or.z4.lt.0d0) then
       znew=zold+dt*k2
     else
       call grid2mk(x1,z4,v,k4)
       znew=zold+(k1+2*k2+2*k3+k4)*dt/6
     endif
   endif

 end subroutine rk4   


!===================================================================
 subroutine rk2bt_substep(xold,zold,dtau,xnew,znew)
! 2nd-order in space/1st order in time Runge-Kutta backward in time
! to find the new position (xnew,znew) of a marker initially 
! at (xold,zold) passively advected by the flow -(u,v)
!===================================================================
   use share; use domain
   implicit none
   real (kind=8), intent(in) :: xold,zold,dtau
   real (kind=8), intent(out) :: xnew,znew
   real (kind=8) :: x1,x2,z1,z2,k1,k2
!===================================================================
   
   x1=xold; z1=zold
   call grid2mk(x1,z1,uextr,k1)
   x2=xold-dtau/2*k1
   call grid2mk(x2,z1,uextr,k2)
   xnew=xold-dtau*k2

   x1=xold; z1=zold
   call grid2mk(x1,z1,vextr,k1)
   z2=zold-dtau/2*k1
   call grid2mk(x1,z2,vextr,k2)
   znew=zold-dtau*k2

 end subroutine rk2bt_substep

!===================================================================
 subroutine rk2bt(i,j,xnew,znew)
! 2nd-order in space/1st order in time Runge-Kutta backward in time
! to find the new position (xnew,znew) of a marker initially 
! at (xold,zold) passively advected by the flow -(u,v)
!===================================================================
   use share; use domain
   implicit none
   integer, intent(in) :: i,j
   real (kind=8) :: xold,zold
   real (kind=8), intent(out) :: xnew,znew
   real (kind=8) :: x1,x2,z1,z2,k1,k2
!===================================================================
   
   xold=xgrid(i); zold=zgrid(j)
   x1=xold; z1=zold

   k1=-uextr(i,j)
   x2=xold+dt/2*k1
   call grid2mk(x2,z1,-uextr,k2)
   xnew=xold+dt*k2

   k1=-vextr(i,j)
   z2=zold+dt/2*k1
   call grid2mk(x1,z2,-vextr,k2)
   znew=zold+dt*k2

 end subroutine rk2bt   


!===================================================================
 subroutine rk4bt_substep(xold,zold,dtau,xnew,znew)
! 4th-order in space/1st order in time Runge-Kutta backward in time
! to find the new position (xnew,znew) of a marker initially 
! at (xold,zold) passively advected by the extrapolated flow (uextr,vextr)
!===================================================================
   use share; use domain
   implicit none
   real (kind=8), intent(in) :: xold,zold,dtau
   real (kind=8), intent(out) :: xnew,znew
   real (kind=8) :: x1,x2,x3,x4,z1,z2,z3,z4,k1,k2,k3,k4
!===================================================================
   
   x1=xold; z1=zold
   call grid2mk(x1,z1,uextr,k1)
   x2=xold-dtau/2d0*k1
   call grid2mk(x2,z1,uextr,k2)
   x3=xold-dtau/2d0*k2
   call grid2mk(x3,z1,uextr,k3)
   x4=xold-dtau*k3
   call grid2mk(x4,z1,uextr,k4)
   xnew=xold-(k1+2d0*k2+2d0*k3+k4)*dtau/6d0

   x1=xold; z1=zold
   call grid2mk(x1,z1,vextr,k1)
   z2=zold-dtau/2d0*k1
   call grid2mk(x1,z2,vextr,k2)
   z3=zold-dtau/2d0*k2
   call grid2mk(x1,z3,vextr,k3)
   z4=zold-dtau*k3
   call grid2mk(x1,z4,vextr,k4)
   znew=zold-(k1+2d0*k2+2d0*k3+k4)*dtau/6d0

 end subroutine rk4bt_substep


!===================================================================
 subroutine rk4bt(i,j,xnew,znew)
! 4th-order in space/1st order in time Runge-Kutta backward in time
! to find the new position (xnew,znew) of a marker initially 
! at (xold,zold) passively advected by the extrapolated flow -(uextr,vextr)
! If one of the four computation points needed lies outside 
! the domain, the new position is found by a 2nd order Runge-Kutta
!===================================================================
   use share; use domain
   implicit none
   integer, intent(in) :: i,j
   real (kind=8) :: xold,zold
   real (kind=8), intent(out) :: xnew,znew
   real (kind=8) :: x1,x2,x3,x4,z1,z2,z3,z4,k1,k2,k3,k4
!===================================================================
   
   xold=xgrid(i); zold=zgrid(j)
   x1=xold; z1=zold

   k1=uextr(i,j)
   x2=xold-dt/2*k1
   call grid2mk(x2,z1,uextr,k2)
   x3=xold-dt/2*k2
   if(x3.gt.dimx.or.x3.lt.0d0) then
     xnew=xold-dt*k2
   else
     call grid2mk(x3,z1,uextr,k3)
     x4=xold-dt*k3
     if(x4.gt.dimx.or.z4.lt.0d0) then
       xnew=xold-dt*k2
     else
       call grid2mk(x4,z1,uextr,k4)
       xnew=xold-(k1+2d0*k2+2d0*k3+k4)*dt/6d0
     endif  
   endif

   k1=vextr(i,j)
   z2=zold-dt/2*k1
   call grid2mk(x1,z2,vextr,k2)
   z3=zold-dt/2*k2
   if(z3.gt.dimz.or.z3.lt.0d0) then
     znew=zold-dt*k2
   else
     call grid2mk(x1,z3,vextr,k3)
     z4=zold-dt*k3
     if(z4.gt.dimz.or.z4.lt.0d0) then
       znew=zold-dt*k2
     else
       call grid2mk(x1,z4,vextr,k4)
       znew=zold-(k1+2d0*k2+2d0*k3+k4)*dt/6d0
     endif
   endif

 end subroutine rk4bt
!========================

 end module markeraux
!########################

